GE ~ power log model from Ettema & Loras, 2009

library(tidyverse)
library(geomtextpath)
library(ggtext)
library(colorspace)
library(sysfonts)
library(showtext)

## use font & icons in ggplots
sysfonts::font_add_google(
    name = "Merriweather Sans", 
    family = "Merriweather Sans"
)
sysfonts::font_add(
    family = "fa-brands", 
    regular = r"(C:\Users\Jem\AppData\Local\Microsoft\Windows\Fonts\Font Awesome 7 Brands-Regular-400.otf)"
)
showtext::showtext_opts(dpi = 600)
showtext::showtext_auto()

## custom theme
theme_set(
    theme_bw(base_size = 12, base_family = "Merriweather Sans") +
        theme(
            plot.title = ggtext::element_textbox_simple(
                size = rel(1.2), lineheight = 1.1
            ),
            plot.subtitle = ggtext::element_textbox(
                colour = "grey35", face = "italic", hjust = 0, 
                lineheight = 1.1, margin = margin(t = 4, b = 4)
            ),
            plot.caption = ggtext::element_textbox(
                colour = "grey35", halign = 1
            ),
            panel.border = element_rect(linewidth = 0.8),
            axis.title = element_text(face = "bold"),
            panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(),
            legend.position = "none"
        )
)

social_caption <- "<span style='font-family:fa-brands'>&#xf09b; &#xe671;</span> jemarnold"

Post 1

Typical range for gross metabolic efficiency (GE) = 15-25%

GE logarithmically increases with power output, so higher fitness = higher power = higher GE

Most of the variance here comes from differences in cadence. The remaining variance comes from random cross-sectional individual differences

## digitized from Ettema & Loras, 2009. Efficiency in cycling: a review Figure 2.d
GE_data <- read.csv(
    r"(C:\R-Projects\data-vis-threads\2026\01-jan\data\ettema-loras-efficiency.csv)"
) |> 
    mutate(
        power = round(power, 1),
        GE = round(GE, 3)
    )

model = lm(GE ~ log(power), data = GE_data)

newdata <- data.frame(
    power = floor(min(GE_data$power)):ceiling(max(GE_data$power))
)

pred <- predict(model, newdata = newdata, interval = "prediction", level = 0.90)

model_data <- cbind(newdata, pred) |>
    tibble() |>
    rename(GE = fit, conf.low = lwr, conf.high = upr)

ggplot(GE_data, aes(x = power, y = GE)) +
    labs(
        title = "Gross Efficiency (GE) increases **logarithmically** with power output",
        subtitle = "Reproduced from Ettema & Lorås, 2009. Fig. 2d.",
        caption = str_glue(
            "**Data**: DOI: 10.1007/s00421-009-1008-7 | **Visuals**: {social_caption}"
        )
    ) +
    theme(
        plot.subtitle = ggtext::element_textbox(margin = margin(t = 4, b = 4))
    ) +
    scale_x_continuous(
        name = "Power Output (W)",
        expand = expansion(mult = 0.01)
    ) +
    scale_y_continuous(
        name = "Gross Efficiency (%)",
        expand = expansion(mult = 0.01)
    ) +
    geom_point(fill = "white", size = 2.4, shape = 21, stroke = 0.5) +
    geom_point(colour = "white", size = 1.6) +
    geom_ribbon(
        data = model_data,
        aes(ymin = conf.low, ymax = conf.high, colour = "GE", fill = "GE"),
        alpha = 0.2
    ) +
    geom_line(data = model_data, aes(colour = "GE")) +
    geom_errorbar(
        data = filter(model_data, power == 250),
        aes(ymin = conf.low, ymax = conf.high),
        linewidth = 0.8, width = 20
    ) +
    annotate(
        "curve",
        x = 310, y = 15, xend = 250, yend = 18,
        curvature = -0.5, arrow = arrow(length = unit(5, "mm"))
    ) +
    annotate(
        "label", x = 350, y = 15, 
        label = "90% prediction intervals\n± 2.5% GE"
    ) +
    geom_point(data = filter(model_data, power == 250), size = 3)